Travel Time Matrix to Maps

Notebook Purpose

This notebook serves to summarize the entire visualization process going from the travel time matrix to the visualizations. That involves the following sections:

1. Travel Time Matrix Wrangling
2. Score Computation
3. Isochrone Computation
4. Dataset Wrangling Part II (NA Insertion)
5. Interactive Visualization
6. Map HTML Exports

0) Useful Libraries

# wrangling/convenience
library(tidyverse)
library(glue)
library(stringr)
library(sf)
library(data.table)

# visualization
library(leaflet)
library(mapview); mapviewOptions(platform = 'leafgl')
#library(ggplot2)
#library(RColorBrewer)
#library(scales)
#library(lattice)

# For pretty knitting
library(lemon)
knit_print.data.frame <- lemon_print
knit_print.tbl <- lemon_print
knit_print.summary <- lemon_print

1) Data Wrangling

Import all dissemination block data

# import dissemination blocks and keep id and pop columns
origins <- fread(file.path("../data/clean", "vancouver_db.csv"))[, .(id, pop)]
origins$pop <- str_replace_all(origins$pop, ',', '')

# change col types
origins$pop <- as.numeric(origins$pop)  
origins$id <- as.factor(origins$id)  

n_origins <- nrow(origins)
paste('Origin Rows: ', n_origins)

[1] "Origin Rows:  15197"

# Peek
head(origins)

Import all amenity data

# import amenities (Cultural/Art facilities)
destinations <- fread(file.path("../data/clean", "vancouver_facilities_2.csv"))

# see summary counts of each amenity
destinations %>% group_by(type) %>% summarise(count = n()) %>% arrange(desc(count))

# clean amenities / filter types to keep 4 most frequent amenities
target_amenities <- c('gallery', 'museum', 'library or archives', 'theatre/performance and concert hall')
destinations <- destinations %>% filter(type %in% target_amenities)
# keep only id and type columns
destinations <- destinations[ , .(id, type)]
# change col types
destinations$type <- as.factor(destinations$type)
destinations$id <- as.factor(destinations$id)  

n_amenities <- nrow(destinations)
paste('Destinations: ', n_amenities)

[1] "Destinations:  354"

head(destinations)

Import the travel time matrix

# import travel time matrix
ttm <- fread(file.path('../data/clean', 'ttm.csv'))

# convert Ids to  factor
ttm$fromId <- as.factor(ttm$fromId)
ttm$toId <- as.factor(ttm$toId)

## Replace travel times less than 1 minute to 1 minute
# This is done to prevent infinity values in the scoring since
# 1 minute is still a reasonable time for trips in the 0 - 1 min range.
ttm$avg_time <- pmax(ttm$avg_time, 1)

# add column for isochrone groupings (between function is inclusive <=)
ttm$time_group <-  with(ttm, ifelse(avg_time < 15, "15",
                          ifelse(avg_time < 30, "30",
                          ifelse(avg_time < 45, "45",
                          ifelse(avg_time < 60, "60",
                          ifelse(avg_time < 75, "75",
                          ifelse(avg_time < 90, "90", "90+")))))))

# how many origins actually have transit accessibility
paste('Origins considered:', round((length(unique(ttm$fromId))/n_origins)*100, 2), '%')

[1] "Origins considered: 94.38 %"

paste('Destinations considered:', round(length(unique(ttm$toId))/n_amenities*100, 2), '%')

[1] "Destinations considered: 97.46 %"

paste('Rows = ', nrow(ttm))

[1] "Rows =  4172482"

# peek
head(ttm)

Import amenity weights but they still have issues:

There are duplicates on IDs 4180 and 4181.

Additionally, of all amenity IDs in the travel time matrix, not all are found in the weight set.

# import amenity weights
amenity_wts <- read.csv('../data/amenity_weights/amenity_wts.csv')

# clean weights
amenity_wts <- amenity_wts[, c('id', 'Index')]
names(amenity_wts) <- c('id', 'weight')
amenity_wts$id <-  as.factor(amenity_wts$id)

# see weight distribution
#plot(density(amenity_wts$weight), main = 'Amenity Popularity Distribution')

amenity_wts %>% group_by(id) %>% summarize(n = n()) %>% arrange(desc(n))

# are all the ttm amenity IDs in the weighted IDs set?
paste('Are all the ttm amenity IDs in the weighted IDs set? (needs to be true for the join to work)')

[1] "Are all the ttm amenity IDs in the weighted IDs set? (needs to be true for the join to work)"

all(unique(ttm$toId) %in% unique(amenity_wts$id))

[1] FALSE

# join when weight index is corrected
#destinations <- left_join(destinations, amenity_wts, by = c('id' = 'id'))

Isochrone Frame and Types

# only consider the nearest amenity
isochrone_frame <- ttm %>%
  group_by(fromId, type) %>%
  summarise(avg_time = min(avg_time))

`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.

# isochrone frame
isochrone_frame$time_groups <-  with(isochrone_frame,
                                   ifelse(avg_time < 15, "15",
                                   ifelse(avg_time < 30, "30",
                                   ifelse(avg_time < 45, "45",
                                   ifelse(avg_time < 60, "60",
                                   ifelse(avg_time < 75, "75",
                                   ifelse(avg_time < 90, "90", "90+")))))))

isochrone_frame <- isochrone_frame[, c(1,2,4)]
head(isochrone_frame)

Import the dissemination block shape file

canada_shape <- st_read("../data/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp", stringsAsFactors = FALSE)

Reading layer `DB_Van_CMA' from data source `C:\Users\Luka\Documents\GitHub\Capstone-Project\data\census2016_DBS_shp\DB_Van_CMA\DB_Van_CMA.shp' using driver `ESRI Shapefile'
Simple feature collection with 15197 features and 27 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 4001643 ymin: 1957237 xmax: 4068119 ymax: 2032663
Projected CRS: PCS_Lambert_Conformal_Conic

# select a greater metropolitan area
metropolitan_area <- "Vancouver"

# filter columns and rows
vancouver_shape <- data.frame(canada_shape[which(canada_shape$CMANAME == metropolitan_area), c(1, 28)])

# id to factor
vancouver_shape$DBUID <- as.factor(vancouver_shape$DBUID)

paste('Rows = ', nrow(vancouver_shape))

[1] "Rows =  15197"

head(vancouver_shape)

2) Score Computation

Notes:

• We don’t scale the data before score computation because we care about the scale of time and standard deviation. We should use their values as is, otherwise scaling will give them equal weighing again when this is a bad assumption.
• Log normalizing the score isn’t important as the score visualization depends on the quantiles taken from the distribution of scores. Since log only shifts values, the visualizations will be identical.

# import score functions
source('Score Modeling/Updated_Score_Functions.R')

score_list <- list()
i <- 1
for (n in c(1, 2, 3, 4)) {
  
  if (n == 4) { n <- NULL }
  score <- sum_score_fxn(ttm, nearest_n = n, weight = FALSE, log_normalize_score = FALSE)
  score_list[[i]] <- score
  i <- i+1
}

`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId', 'type'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId', 'type'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'fromId'. You can override using the `.groups` argument.

scores_long <- rbindlist(score_list) %>% arrange(fromId, nearest_n)

scores_long

NA

3) Dataset Wrangling Part II (NA Insertion)

Each origin(fromId) should have x different scores based on the possible equation below:

In the current file 1 weight options * 4 amenity options * 4 nearest options = 16

In the future 2 weight options * 4 amenity options * 4 nearest options = 32

Many zones however don’t have any routes, such as provincial parks which aren’t linked to via transit. To visualize them in the data we need to attribute these cases an NA value. Since these rows are missing we’ll use the following code to re-add those NA values.

# target zone count
x <- 16

# check for the inconsistency in values
counted <- scores_long %>% group_by(fromId) %>% summarize(n = n()) %>% arrange(n) 

# check for current unique counts to make sure it's multiples of the type number
paste('Unique counts: '); unique(counted$n)

[1] "Unique counts: "
[1]  4  8 16

# get expected rows by multiplying unique IDs by x
n_fromIds <- uniqueN(scores_long$fromId)
N <- n_fromIds*x

paste(glue('{nrow(scores_long)} of {N} rows filled ({round((nrow(scores_long) / N)*100, 2)}%)'))

[1] "229128 of 229488 rows filled (99.84%)"

paste(N - nrow(scores_long), 'to fill.')

[1] "360 to fill."

# function for NA grid expansion 
NA_grid_maker <- function(id, df, isochrone = FALSE) {
  
  all_amenities <- as.character(unique(df$type))
  
  # get missing amenities by indexing the fromId and keeping only unique types
  missing_amenities <- setdiff(all_amenities, unique(df$type[df$fromId == id]))
  
  if (isochrone == FALSE) {
    # create NA rows to append via expand.grid (creates a row for every factor combination)
    NA_rows <- expand.grid('fromId' = id,
                         'type' = missing_amenities,
                         'weight' = as.character(unique(df$weight)),
                         'nearest_n' = as.character(unique(df$nearest_n)),
                         'score' = NA, 
                         stringsAsFactors = TRUE)
  } else {
    # create NA rows to append via expand.grid (creates a row for every factor combination)
    NA_rows <- expand.grid('fromId' = id,
                         'type' = missing_amenities,
                         'time_groups' = NA, 
                         stringsAsFactors = TRUE)
  }
  
  NA_rows
}

# function for filling table with NA values
NA_table_filler <- function(df, custom_idx = NULL, isochrone = FALSE) {
  
  # count each fromId occurence
  fromId_counts <- df %>% group_by(fromId) %>% mutate(n = n())
  
  if (is.null(custom_idx)) {
    # create a fromId array using Ids that don't meet the [x] count requirement
    id_arr <- array(unique(fromId_counts[fromId_counts$n < x, ]$fromId))
  } else {
    id_arr <- custom_idx
  }
  
  # get rows
  filler_rows <- rbindlist(apply(id_arr, MARGIN = 1, FUN = NA_grid_maker, df = df, isochrone = isochrone))

  # append and order
  df <- rbindlist(list(df, filler_rows), use.names = TRUE) 
  
  if (isochrone == FALSE) {
    df <- df %>% arrange(fromId, type, nearest_n, weight)
  } else {
    df <- df %>% arrange(fromId, type)
  }
  
  df
}


partial_filled_scores_long <- NA_table_filler(scores_long)

Although, recall that not all origins were used in the score computation, so we need to add those NA values as well.

Now lets add population data to the scores and isochrone frame

## Export checkpoint
write.csv(filled_scores_long, 'June2_scores_long.csv', row.names = FALSE)
write.csv(filled_isochrone_frame, 'June2_isochrone_frame.csv', row.names = FALSE)

4) Interactive Visualization

# join factor and geometry data 
scores_viz_frame <- left_join(vancouver_shape, filled_scores_long, by = c('DBUID' = 'fromId'))
isochrone_viz_frame <- left_join(vancouver_shape, filled_isochrone_frame, by = c('DBUID' = 'fromId'))

# convert back to sf object
scores_viz_frame_sf <- st_as_sf(scores_viz_frame)
scores_viz_frame_st <- st_transform(scores_viz_frame_sf, crs = 4326)

isochrone_viz_frame_sf <- st_as_sf(isochrone_viz_frame)
isochrone_viz_frame_st <- st_transform(isochrone_viz_frame_sf, crs = 4326)

sort(unique(isochrone_viz_frame_st$time_groups))

[1] "15"  "30"  "45"  "60"  "75"  "90"  "90+"

map_maker_isochrone(isochrone_viz_frame_st, 'gallery')

[1] "Current Map:"
Gallery Transit Isochrone

5) Map HTML Exports

for (amenity in unique(scores_viz_frame_st$type)) { 
  
  # 4 isochrone maps
  map_maker_isochrone(data = isochrone_viz_frame_st, amenity)

  for (weight in unique(scores_viz_frame_st$weight)) {
    for (nearest_n in unique(scores_viz_frame_st$nearest_n)) {
      
      # 16 others
      map_maker_scores(data = scores_viz_frame_st, amenity, weight, nearest_n)
      
    }
  }
}